 rm(list = ls())
 setwd("E:/ASM/snp-noWASP-nocountdiff-onlyWGS-T&N")
file=read.table("E:/ASM/countdiff只考虑单发/interested_position/unique_id_61positions_in_M49-X2B_WGS.txt",head=T,sep="\t")

positions=unique(file$unitID)
positions=as.character(positions)
sel=list.files(pattern=".snp")
positions
i=1
file=read.table(sel[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final=final[complete.cases(final[,ncol(final)]),]		#去掉空值
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
final$con_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)
final$case_reads=as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
#sels=which(final$con_reads>=10 | final$case_reads>=10)
sels=which(final$total_reads>=15)
final=final[sels,]
name=names(final)
final=final[,name[-c(8:10)]]
final$DU=abs(final$normal_var_freq-0.5)
final$DA=abs(final$tumor_var_freq-0.5)
cols=ncol(final)
final$p.diff=lapply(1:nrow(final),function(x) fisher.test(rbind(c(final[x,1],final[x,2]),c(final[x,4],final[x,5])))$p.value)
final$p.biasA=lapply(1:nrow(final),  function(x) binom.test(c(final[x,4],final[x,5]),p = 0.5)$p.value)
final$p.biasU=lapply(1:nrow(final),  function(x) binom.test(c(final[x,1],final[x,2]),p = 0.5)$p.value)

colnames(final)=c(paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID",paste(gsub(".snp","",sel[i]),c("DU","DA"),sep="_"),paste(gsub(".snp","",sel[i]),c("p.diff","p.biasA","p.biasU"),sep="_"))
}

i=2
file=read.table(sel[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final1=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1=final1[complete.cases(final1[,ncol(final1)]),]		#去掉空值
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
final1$con_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)
final1$case_reads=as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
#sels=which(final1$con_reads>=10 | final1$case_reads>=10)
sels=which(final1$total_reads>=15)
final1=final1[sels,]
name=names(final1)
final1=final1[,name[-c(8:10)]]
final1$DU=abs(final1$normal_var_freq-0.5)
final1$DA=abs(final1$tumor_var_freq-0.5)
cols=ncol(final1)
final1$p.diff=lapply(1:nrow(final1),function(x) fisher.test(rbind(c(final1[x,1],final1[x,2]),c(final1[x,4],final1[x,5])))$p.value)#计算各种统计值
final1$p.biasA=lapply(1:nrow(final1),  function(x) binom.test(c(final1[x,4],final1[x,5]),p = 0.5)$p.value)
final1$p.biasU=lapply(1:nrow(final1),  function(x) binom.test(c(final1[x,1],final1[x,2]),p = 0.5)$p.value)
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID",paste(gsub(".snp","",sel[i]),c("DU","DA"),sep="_"),paste(gsub(".snp","",sel[i]),c("p.diff","p.biasA","p.biasU"),sep="_"))
final=merge(final,final1,by="unitID",all=T)

}
final2=data.frame(lapply(final,as.character),stringsAsFactors=FALSE)